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Wc study the dynamical correlation functions of the Richardson pairing model (also known as the 
reduced or discrete-state BCS model) in the canonical ensemble. We use the Algebraic Bethe Ansatz 
formalism, which gives exact expressions for the form factors of the most important observables. 
By summing these form factors over a relevant set of states, we obtain very precise estimates 
of the correlation functions, as confirmed by global sum-rules (saturation above 99% in all cases 
considered). Unlike the case of many other Bethe Ansatz solvable theories, simple two-particle 
states are sufficient to achieve such saturations, even in the thermodynamic limit. We provide 
explicit results at half-filling, and discuss their finite-size scaling behavior. 



I. INTRODUCTION 

Possibly the most remarkable incipient property asso- 
ciated to a fermionic gas is its instability to the pairing 
phenomenon. Under an arbitrarily weak attractive force, 
which can originate in such a simple process as coupling 
to phonons, the gas will develop an instabiUty towards 
the formation of Cooper pairs, this process forming the 
basis of the BCS theory of superconductivity^, with fur- 
ther remarkable consequences such as the Meissner and 
Joscphson effects. Within BCS theory, single-particle 
excitations are suppressed by the superconducting gap, 
which is obtained from the solution of a variational ansatz 
for the wavefunction defined in the grand-canonical en- 
semble (electron number is by construction not conserved 
anymore) and in the thermodynamic limit (since a con- 
tinuous energy band is assumed). 

While these mean-field approaches successfully de- 
scribe the experimental features of the traditional bulk 
superconductors, recent experiments have also consid- 
ered metallic nanograins^'^ in which the level spacing is 
finite and of the same order as the superconducting gap, 
and in which some Coulomb blockade effects could oc- 
cur in view of the finite charging energy of the grain. 
These studies open the door to many interesting fur- 
ther questions not answerable within BCS theory, and 
require mesoscopic effects to be encompassed back into 
the model. One way to approach this problem is to use 
the so-called reduced BCS model, defined by the Hamil- 
tonian 

JV N 

HbCS= Y^acrCaa-g Yl <i+cl-C0-C0+, (1) 

c=i a, 13=1 

which was introduced by Richardson in the early 1960's 
in the context of nuclear physics^. The model describes 
(pseudo) spin-1/2 fermions (electrons, nucleons, etc..) 
in a shell of doubly degenerate single particle energy 
levels with energies ea/2, a = 1,...N. Ca^a are the 
fermionic annihilation operators, cr = +, — labels the de- 
generate time reversed states (i.e. spin or isospin) and 



g denotes the effective pairing coupling constant. De- 
spite its simplified character (the interaction couples all 
levels uniformly), the model does have a number of ad- 
vantages as compared to BCS theory. First of all, it can 
be solved within the canonical ensemble (fixed number 
of electrons), a situation which is relevant for isolated 
nanograins. Second, and rather remarkably for an ex- 
actly solvable model, it remains solvable for an arbitrary 
choice of parameters, and can thus provide quantitative 
predictions for various situations obtained by considering 
various choices of the set of energy levels Ca (both their 
number, and their individual value), coupling g, and fill- 
ing. Besides mesoscopic superconductivity, this model 
and its solution also find applications in other fields (see 
the reviews [5,6] for some applications outside of con- 
densed matter physics). 

The nature of the electronic states in a metallic 
nanograin can conceivably be probed in a number of dif- 
ferent experiments. Electronic transport through such a 
grain could be studied by attaching either metallic or su- 
perconducting leads. The observable I-V characteristics 
or Joscphson currents would be theoretically obtainable 
from correlation functions within the grain. Such corre- 
lations, however, are not easily obtainable from the basic 
exact solution of the model, which focuses on wavefunc- 
tions but does not allow to make direct contact with the 
dynamics of observables. The history of the study of 
correlations in the Richardson model is however already 
rather rich. Richardson himself in 1965^ derived a first 
exact expression for static correlation functions. In a 
significant development, Amico and Osterloh^ proposed 
a new method to write down such correlations explic- 
itly but with results limited to system sizes of up to 16 
particles. A major simplification was then proposed by 
Zhou ct al.^'^'^: using the Algebraic Bethe Ansatz (ABA) 
and the Slavnov formula for scalar products of states^ ^, 
they managed to obtain the static correlation functions 
as sums over determinants of Np x Np matrices. These 
expressions have been further simplified by us in a previ- 
ous publication, where they were evaluated numerically 
for a particular choice of the energy levels ea^'^ allowing 
to describe the crossover from mesoscopic to macroscopic 
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physics, going beyond previous results limited to fewer 
particles^'^"^. We mention that a different approach, valid 
in the case of highly degenerate £„, is also available^^. 

Despite all these developments, up to now none of these 
approaches has been adapted and used to calculate dy- 
namical correlation functions, important to quantify the 
response of a physical system to any realistic experimen- 
tal probe. In this paper we use a method that mixes 
integrability and numerics and is similar to that used by 
some of us to study the dynamical correlation functions^^ 
and entanglement entropy-'^^ in spin-chains, as well as dy- 
namical correlations in Bose gascs^^. In short, the ABA- 
CUS method^^ consists in using the exact knowledge of 
the form factors of physical obscrvables as determinants 
of matrices whose entries are the unknown Richardson 
rapidities. For any given state, these rapidities are cal- 
culated by solving the Richardson equations. The cor- 
responding form factors arc then evaluated. Finally we 
need to sum over all these contributions, but this can be 
done by searching for the states in the Hilbert space that 
contribute more significantly to the correlation function. 
This is done by optimizing properly the scanning of the 
Hilbert space^^ and the accuracy of the result is kept 
under control by checking the values of the sum-rules. 
However, we will see that for the Richardson model, this 
scanning is particularly easy since the two-particle states 
dominate the sum even when increasing the number of 
particles (in strong contrast with what happens for other 
models^"'''^^). This property is clearly connected with the 
mean-field character of the model in the thermodynamic 
limit and the subsequent suppression of quantum fluctua- 
tions. However at finite (sufficiently low) number of par- 
ticles the effect of quantum fluctuations can be revealed 
by small, but maybe measurable, multi-particle chan- 
nels. We mention that another variation on the ABA- 
CUS logic has recently been used for the calculation of 
out-of-equilibrium observables in the pairing model sub- 
jected to a quantum quench, but in that case the class of 
excitations contributing was much wider 

The paper is organized as follows. In Sec. II we dis- 
cuss the model and its general properties. In Sec. HI we 
recall the Algebraic Bethe Ansatz approach to the model 
and we proceed to some simplifications of the determi- 
nant expressions. In Sec. IV we introduce the sum rules 
for the considered correlation functions and we present 
the selection rules for the form factors that will help in 
the numerical computation by reducing the number of 
intermediate states we have to sum over. In Sec. V all 
the correlation functions are explicitly calculated at half- 
filling. We report our main conclusions and discuss open 
problems for future investigation in Sec. VI. Appendix A 
reports the technical details on how to solve the Richard- 
son equations for any excited state while Appendix B 
looks at the details of the strong coupling expansion of a 
studied correlation function in order to explain its scaling 
behavior. 



II. THE MODEL 

As it is written in Eq. (1), singly- as well as doubly- 
occupied levels arc allowed. The interaction however cou- 
ples only doubly-occupied levels among themselves, and 
due to the so-called blocking effect^'^, unpaired parti- 
cles completely decouple from the dynamics and behave 
as if they were free. We will denote the total number of 
fermions as Nf, and the total number of pairs as Np. Due 
to level blocking, we will thus only consider Nf = 2Np 
paired particles in N unblocked levels, keeping in mind 
that we could reintroduce blocked levels later if needed 
in the actual phenomenology desired. In terms of pair 
annihilation and creation operators 

ba — Ca — Cq.-(- 6]j — C|!j,_j_C|!j_ , (2) 

the Hamiltonian is 

N N 

H=Y^e^bih^~gY.hlbp, (3) 

a=l a,/3=l 

and Ua = '2b\j}a is the number of particles in level a. 
The pair creation and annihilation operators satisfy 

the commutation relations 

hi] = 5„^(1 - 2bib^) , [6«, hp] = [bi,bl\ = . 

(4) 

The term 2b\jDa in the first commutator makes the model 
different from free bosons and therefore non-trivial. 
Using the pseudo-spin realization of electron pairs 
- hib^ - 1/2, S-=b^, 5+ = hi, the BCS Hamilto- 
nian becomes (up to a constant) 

N N 

H=Y,e^Sl-gY, SlS-p. (5) 

The operators S*^'^ obey the standard su{2) spin algebra 
and so the Hamiltonian (5) describes a spin-1/2 magnet 
with long-range interaction for the XY components in 
a site-dependent longitudinal magnetic field £„. Such 
a magnetic Hamiltonian is known in the literature as a 
Gaudin magnet^^. An important relation is 

A. Grand-canonical BCS wavefunction 

In the grand-canonical (GC) ensemble the ground state 
wavefunction is the BCS variational ansatz 

\GS)=J{{u^ + e^'^-Vo.bi)\Q), ul+vl = l, (7) 

a 

where the variational parameters Ua and u„ are real and 

(pa is a phase which, it turns out, must be a-indepcndent. 
\GS) is not an eigenstate of the particle number operator 
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Nf and the average condition (iV/) = Nf determines the 
GC chemical potential Likewise, the commonly used 
definition 



(8) 



for the superconducting gap makes sense only in a GC 
ensemble, since {be) is zero when evaluated at fixed par- 
ticle number. The variational parameters are obtained 

as 
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1- 



V(ea-M)" + |AGcP 
where /U is the GC chemical potential. 



(9) 



B. Canonical description and Richardson solution 



Throughout the paper we will refer with latin indices to 
the rapidities and with greek ones to the energy levels. 
The total energy of a Bethe state is, up to a constant, 



E. 



{w} 



E 



Wi 



(14) 



For a given A'' and Np the number of solutions of the 
Richardson equations is (jy ) , and coincides with the di- 
mension of the Hilbert space of Np pairs distributed into 
N different levels, i.e. the solutions to the Richardson 
equations give all the eigenstates of the model. 

Note that one could also use the pseudovacuum to be 
the state fully polarized along the z axis (instead of— 5) 
This leads to slight differences in the expressions, which 
one can resolve by comparing references [9] and [10]. 



The exact solution (i.e. the full set of eigenstates and 
eigenvalues) of the Hamiltonian (1) in the canonical en- 
semble was derived by Richardson'*. The model can be 
encompassed into the framework on integrable models^^ 
and is tractable by means of algebraic methods^'^"'^^'^"*. 
We review here only the main points of this solution. 

In the ABA, eigenstates are contructed by applying 
raising operators on a so-called reference state (pseu- 
dovacuum). We here choose the pseudovacuum (in the 
pseudo-spin representation) to be fully polarized along 
the —z axis 



5^|0) = -^|0), Va. 



(10) 



In the pair representation, this state thus corresponds to 
the Fock vacuum. Eigenstates with Np pairs are then 
characterized by Np spectral parameters (rapidities) Wj, 
and take the form of Bethe wavefunctions 



N„ 



\{wj}) = i[c{wkm. 



(11) 



fe=i 



The operators C, together with operators A, B, V defined 
as 

9 ~,'^k-ea ~,Wk- ea 



N 



Si 



a=l 
N 



a=l 



Wk - <^a 



V{wk) 



1 -^^{12) 
z-^ ,,,, _ ^_ \ ' 



9 ^Ti - ^ 



obey the Gaudin algebra, which is the quasi-classical 
limit of the quadratic Yang-Baxter algebra associated to 
the gl{2) invariant i?-matrix (we refer the readers to [10] 
for details). 

The wavefunctions (11) are eigenstates of the trans- 
fer matrix, and thus of the Hamiltonian (1), when the 
parameters wj satisfy the Richardson equations 



1 AT 1 

- = J2— — E 



g '■ — ' w. 



j = l,...,Np. (13) 



III. CORRELATION FUNCTIONS AND 
ALGEBRAIC BETHE ANSATZ 

We are interested in the dynamical correlation func- 
tions of the form 



Go„,(t) = 



{GS\Oi{t)Of,mGS) 
{GS\GS) 



(15) 



at fixed number of pairs Np. Here Oa stands for a 'local' 

operator in the Heiscnberg picture (i.e. depending on a 
single energy level a). We consider equal to or 

By inserting the complete set of states with {w} 

a set of rapidities solution to the Richardson equa- 
tions, and using the time evolution of the eigenstate, we 
can rewrite the dynamical correlation function as the sum 



E 

{w} 



{{w}\Oo,\GSr{{w}\Op\GS)e'^-' 
{GS\GS){{w}\{w}) 



(16) 



where form factors and norms are obtainable by Alge- 
braic Bethe Ansatz^"'^'^ and will be discussed in the next 
subsection. The frequency ojw is just tOyj = — Eqs — 
n{Myj — Np), where the energies E/^ and Eqs are given in 
equation (14). /i is the chemical potential (needed only 
for ) and M^, is the number of rapidities of the state 



More easily, we can write i 



E,, 



■ Eo,w, whci 



Eq^w is the lowest energy state with the same number of 
rapidities as Notice that for correlation func- 

tion, we only need states with = Np, and so the 
chemical potential term is absent while for {S~^S~), we 
need only states with Np — 1 rapidities. 

The most relevant physical observablcs are clearly 
global ones, when the sum over the internal energy levels 
is performed. We will consider diagonal and global cor- 
relation functions, whose static counterparts for are 
the diagonal^^ and off-diagonal order parameters. We 



4 



consider in the following the three global correlators 



N 



a=l 
N 



Gi-(t) = E 



a=l 
N 



GlUt) = E 



{GS\SUt)Sm\GS) 
{GS\GS) 

{GS\S+it)S-mGS) 
{GS\GS) 

{GS\S+it)SUO)\GS) 



{GS\GS) 



(17) 
(18) 
(19) 



We will mainly study them in frequency space, since their 
structure is less complicated, and we will only plot a few 
examples in real time. 

We stress that while the level-resolved correlators de- 
pend strongly on the choice of the energy levels e^, for 
the global ones it is expected that most of the qualitative 
features and several quantitative ones are not affected by 
the choice of the model. Thus our results, even if ob- 
tained for a specific choice of , should display the main 
features of the dynamical correlation functions for a wide 
variety of Richardson models. 



from which the norms of states simply follow from v ^ w 
as I III'} IP = detAT^ G with a Gaudin matrix 



N 



G, 



ah 



{Va - VbY 



a ^ 6, 
(23) 

recovering Richardson's old result^. 

The key point is that any form factor of a local spin op- 
erator between two Bethe eigenstates can be represented 
via (12) as a scalar product with one set, e.g. {v} not 
satisfying the Bethe equations, for which Slavnov's for- 
mula is applicable. This has been explicitly worked out 
in Ref. [9]. For containing respectively Np + 1 

and Np elements, the non-zero form factors are: 



{{w)\S-\{v}) = {{v}\Sl\{w}) 

detAr^+ir(Q;, {w},{w}) 



nf=lK - £„) n6>aK - Wa) n6<a(«b " ^a) 

and, for both {w} and {v} containing Np rapidities 



(24) 



A. Algebraic Bethe Ansatz and Form Factors 

The starting point to calculate correlation functions 
with the Algebraic Bethe Ansatz is having a representa- 
tion for the scalar products of two generic states defined 
by Np rapidities {Np Cooper pairs) 



{{w}\{v}) = {0\]jB{wb)]lC{vam, (20) 

6=1 a=l 

when at least one set of parameters (e.g. wi, but not 
Va) is a solution to the Richardson equations. Following 
standard notations, C is the conjugate of the operator B. 
Such a representation exists, and is known as the Slavnov 
formula^ ^, which for the case at hand specifically reads^ 



(MIW) = 



Ub<a(wb - Wa) UaKbi'^b " Va) 
XdetN,J{{Va},{Wb}), (21) 



where the matrix elements of J are given by 



Jab = 



Vb - Wb 



N 

E 



^« ~ '^b V^l ~ ea)(Wb - £«) 

T^^T^^h (22) 



[Va - Vc){Wb -Vc) j ' 



a=l " 

detjv, (|7;(M, H) - Q{a, {w}, {v})) 



Ub>aiWb-Wa)Ub<ai^b-Va) 

with the matrix elements of T given by {b < Np + 1) 

Np+l / N ^ 

Tabia)= U(-e-Vb)(^^^^—^^^^^ 



(25) 



{Vb - Wc){Wa - Wc 



TaNp + lia) 



{Wa - ea)2 



, Qab(a) 



ric^bC^c - ^ 6) 
{Wa - Ca) 



1 ■ 



Above, Tz is the Np x Np matrix obtained from T by 
deleting the last row and column and replacing A'p -|- 1 
by Np in the matrix elements. Here it is assumed that 
both {wa} and are solutions to Richardson's Bethe 
equations. However, the results are still valid for if 
only {wb\ satisfy the Bethe equations. 

When approaching a bifurcation point (see App. A) 
in the solutions of the Richardson equations, some indi- 
vidual terms in the sum defining the matrix elements of 
T tend to diverge (because Wb — > Eq for some h and a). 
Those divergences cancel out when the sum is taken, but 
they can still lead to large numerical inaccuracies. By 
using the Richardson equations, it is possible to elimi- 
nate such potentially problematic terms and rewrite the 
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matrix T as 



Wc 



Vb) 



Tab = 



Vb 



E 



{Vb - Vc) 



-E 



{Vb - Wc) 



(26) 

where potentially diverging terms have been removed. 
We stress once again that such a formula holds only if 
both {wa} and {va} are solutions to the Richardson equa- 
tions with the same g. This expression has been obtained 
before in Ref. 19. 



IV. SUM RULES AND SELECTION RULES 

The standard way to assess the accuracy of a numer- 
ical calculation in which we discard part of the states 
consists in using sum rules, e.g. summing of all the con- 
tributions independently of the energy of the intermedi- 
ate state. When summing over all states, wc always get 
static quantities, and, in the present case, they can be 
obtained by very simple considerations. 

For G^, we have the sum rule 



N 



EE 

a=l {v} 



\{{v}\S^^\GS)\^ 
{GS\GS){{v}\{v}) 



N 

E 



{GS\{S^J^\GS) 



{GS\GS) 



N 
1 



(27) 



Only the (^^ ) states with the same number of rapidities 
as the ground state contribute to this correlation func- 
tion. Instead for Gi^_ we have 



N 

E 



{GS\S+S-\GS) N_ 
{GS\GS) ~ 2" 



N 



(28) 



as easily shown by using Eq. (6). Here, only the 

states with one less rapidity than the ground state con- 
tribute to this correlation and to the similar one contain- 
ing off-diagonal terms G°^_ . In this last case we have 



N 

E 



{GS\StS-p\GS) 
{GS\GS) 



od •) 



(29) 



that is the off-diagonal order parameter, which can be 
obtained by the solution of the Richardson equations 
for the ground state and using the Hellmann-Feynman 
theorem^^. 

We will see in the following that the two-particle states 

will give most of the contribution to the correlation func- 
tions, always saturating the sum-rules to more than 99% 
accuracy. We show in the following subsections that 
some selection rules imply that only two-particle states 
have non-zero contribution to the correlation functions 



for .g = and for g ^ oo. Although at intermediate cou- 
plings this set of states docs not give 100% saturation of 
the sum rules, these two limits clearly give insight as to 
why they remain extremely dominant in every regime. 



A. Weak coupling regime 

In the non-interacting g = Q limit, the fixed Np eigen- 
states are quite naturally described by placing the Np 
Cooper pairs (flipped pseudo-spins) in any of the i^j^ ) 
possible sets of Np energy levels picked from the N avail- 
able ones. This translates into a representation in terms 
of rapidities given by setting the Np rapidities to be 
strictly equal to the energies Cq, of the Np levels occu- 
pied by a pair. Since 



lim C{u) — lim 



St 



(30) 



the states built in such way will have diverging norms and 
form factors, but it remains possible to describe the limit 
correctly because in the ratio of form factors and norm 
the two divergences cancel. Since any of these states 
is an eigenvector of every operators, at 5 = the 
only contributions to the S'^ correlations come from the 
ground state to ground state form factor {GS\ \GS). 
In a perturbative expansion^^ in g, it is easy to see that at 
first order, the corrections to the ground state comes only 



from states 



{w} = |ec('j 1^ differing from it by at 

most one rapidity (in the 5 — >■ limit). They constitute 
the full set of two-particle states, obtained by creating 
a "hole" and a "particle", i.e. moving a single Cooper 
pair (rapidity) in the ground state to any of the available 

unoccupied states. 

At 5 = 0, the {{v}\S^ \GS) form factors are non- 



zero whenever 



{v} = |e„i...ea„^_i|^ is obtained by re- 
moving a single rapidity from the Np pairs ground state 
\GS = {wi = ei, ... uiTVp = ^Np})- These states can also 
all be thought of as two-particle states in the Np — 1 pairs 
sector, since they can all be generated by moving a single 
rapidity in the Np — 1 ground state. 

In the specific case of half filling {Np = N/2), treat- 
ing every two-particle excitation means that only N'^/A 
states are needed, out of the full (^^3) dimensional 
Hilbert space. Quite naturally, when a non-zero coupling 
is included these states might not be sufficient anymore, 
but as will be shown in the next two sections, for g ^ 00, 
only this set of states is once again needed to compute 
every non-zero form factor of local spin operators. 



B. Strong coupling {g 00) regime 

Yuzbashyan et al.^'^'^* showed that the solutions to the 

Richardson equations arc such that in the g ^ 00 limit, 
a number Nr of the rapidities will diverge as Wi ~ Gig + 
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0{g^). The coefficients Ci are given by the roots of 
the appropriate Laguerre polynomiaP^ 



-1-N-2{N,.-Np) 



(Ci) = 0. 



(31) 



In the infinitely large coupling limit, the impact of the 
diverging rapidities is well defined. In fact, lim C{u) oc 

u—>oo 

= is the total spin raising operator. In this 

a 

limit, the coupling term g''^S^S^ = gSfg^S^g^ com- 

afi 

pletely dominates the Hamiltonian and so its eigenstates 
become eigenstates of the Sj^^j operator too. Through 
simple energetic consideration^^'^'', one can then easily 
show that the number of diverging rapidities is related 
to the eigenvalue J( J + 1) of Sf^j by the relation 



N 



(32) 



A state defined by Np rapidities, of which are infinite, 
therefore belongs to a subspace of the total Fock space 
defined by eigenvalues of S^^^ and S^^j given by quantum 
numbers m = -N/2 + Np and J = N,. - N,, + N/2. This 
allows us to derive explicit strong-coupling selection rules 
for the various form factors used in this work. 

Since a given eigenstate, in this limit, is built out of a 
linear superposition of the various spin states with fixed 
J and m, it can be decomposed onto the joint eigenbasis 
of spin a combined with the various multiplets emerging 
from the addition of the remaining N — 1 spins ^ . 

For a given (fixed degeneracy index k) N — 1 spins 
multiplet with magnitude J', one can write the highest 
weight state as 

k,J = J' + ^,m=J' + ^^ = Ita) \k, J', m' = J') . 

(33) 

The n-times repeated action of Sf^t = +<S'jy_j on this 
state will generate the eigenstates 

k,J'+^,m'+^-n^ = C^\W®\k,J',J' -n) 

+ \k, J', J'-n + l), 

where we do not need to explicitly specify the Clebsh- 
Gordan coefficients. Combining this multiplet with a 
spin i also gives rise to a second set of states given by 
J = J' — 1/2 and the corresponding allowed values of 
m G {— J, — J + 1, ... J — 1, J}. These states are easily 
constructed by making them orthogonal to the previously 
found ones, i.e. 



k,J'-^,m' + ^-n^= -CJ Ita) ® \k, J', J' - n) 
+ Cnia)<^|fc, J', J'-n+l). (34) 



Any general state with fixed J and m can therefore 
have contributions coming from J+l/2orJ— 1/2 mul- 
tiplets of the N — 1 excluded spins, i.e.: 



^Ak\k,J,m) = Y,(Bi\W 

k k ^ 

+B2 \ia) 



k,J~\,m-\ 



k,J- 



k,J+ \,m- ^ 



1 



k,J + 2'"^+ 2 



.(35) 



Given this form, it is trivial to see that the application 
of or S~ on any of these states will result in 



k,J-\,m-\ 



—B2 \ia) <: 

+Bl Ita) C 
-B4 14a) ' 



k,J-^,m+^ 
k,J + \,m- ^ 
k,J + 



m) oc ^ 



Bi \ia) 



+Bl \ia) ' 



k,J-\,m-\ 



k,J+ \,rn-^ 



Consequently, the form factors for 

^k' {k', J", m"| j 55 Ak \k, J, m) j , (36) 

can exclusively be non-zero if J" € {J — 1, J, J + 1} and 
m" = m. 

Similarly, the S~ form factor 

^ A',, {k', J", m"| j S- ^ Ak \k, J, m) j , (37) 

is non-zero only if J" G {J — 1, J, J + 1} and m" = m—1. 

Since, as we pointed out earlier, the value of J is related 
to the number of diverging rapidities, these selection rules 
translate into selection rules for the total number of ra- 
pidities and the number of diverging rapidities. 

The 5*^ form factor is non-zero for m" = m and 
therefore for a total number of rapidities in the in- 
termediate states {Np) given Np = Np. Using this 
fact, the selection rule on J" and Eq. (32), we easily 
find that the only contributions come from states with 

N^' = {Nr-l,Nr,Nr + l}. 

Similarly, for S~ , we find that A''^' = Np — 1 and the 
number of diverging rapidities must be given by A^" = 

{Nr,Nr -l,Nr- 2}. 
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For the specific case of ground state (iV,. = Np diverg- 
ing rapidities) expectation values, non-zero contributions 
to 5^ correlations are found exclusively for intermediate 
states with either Np or Np — 1 diverging rapidities since 
TV" = Np -I- 1 > Np is impossible. Identically, for S~ , the 
only possible cases are given by an intermediate state 
with Np — 1 or Np — 2 diverging rapidities. 

Finally, for the form factors involving only the 
ground state, one can use the fact that: 



lim \GS) ^ lim \\ NpC{v,) |0) cx (5, 



|0) 



I N-l \ 

{il, ... QNp-l} 
( Np ) 

{ai, ... ajvp} 



(38) 



where |{tai ■•■ Tqm}) is simply the sum over 

all possible states containing M up spins picked out of 
the iV — 1 levels which exclude level a. The action of 
on this state is trivially given by multiplying by 1/2 
while adding a minus sign to the second sum, which in 
the end allows us to simply prove that: 



lim {GS\ 5^ \GS) cx - 

g->-oo 2 

1 (A^-1)! 



N-l 
N„ 



2 {N - Np - l)\{Np - 1)1 





rN- 










2Np- 


N 





(N - Np)Np 



(39) 



One immediately sees that the half-filling 2Np — N 
case has the peculiar feature of having these 'ground 
state to ground state' form factors go down to zero in 
the strong limit coupling. 



C. Correspondence between g — and g ^ oo 
eigenstates 

1. General algorithm 

An algorithm relating the 5 = structure of a state to 
the number of diverging rapidities it will have at 5 cx) 
was already proposed in Ref.'^°^'^^. It necessitates the 
evaluation of various quantities for every possible parti- 
tions of the N levels into 3 disjoint contiguous sets of 
levels. An equivalent result can be obtained through the 
simple following algorithm, which was discussed before 
in Ref. 20. 

By splitting any 5 = configuration of rapidities into 
blocks of contiguous occupied and empty states, one can 
simply obtain the number of diverging roots. Fig. 1 
shows some examples of this construction. Every circle 



Pi = 3 



Hi = 4 



P3 



o m Hs ^ 1 

1 1 



Pi = 5 



o 
o 
o 

P2=2 J 



Hi = 3 



Al =2; Nr=4 



Pi = 1 



m 



P2 = 2 



Pa = 2 



P4 = 2 



Hi = 2 



H2 = 2 



Ha = 2 



o |o] H4 = 1 

Al =0; A2=0; 
A3 = D; Nr = 1 



FIG. 1: Construction of the contiguous blocks necessary to 
establish the <? = 0, </ — >■ cxa correspondence. 



represents a single energy level and the blackened ones 
are occupied by a Cooper pair at 5 = 0. 

We then label the various blocks according to the fol- 
lowing prescription. The highest block of rapidities is 
labelled by index i = 1 and contains Pi rapidities. The 
block of unoccupied states right below it will be also la- 
belled by i = 1 and contains Hi empty levels. We con- 
tinue this labeling by defining i-2(^^2) as the number of 
rapidities (unoccupied states) in the next block until ev- 
ery single one of the Ni, blocks has been labelled. In the 
event that the lowest block (i = N^) is a block of rapidi- 
ties, as is the case in the middle example in Fig. 1, we 
set Hmi, = 0. 

The number of diverging rapidities is then simply given 
by: 

Nr = [Pn, + ^w.-i - MHPn, + ^iv,-i, i^ivJ] (40) 
with the Ai terms defined recursively as: 

A, = [P, + A,_i - Min(P, + A,^i,Hi)] (41) 
with = 0. 

This can be thought of as a 'dynamical' process which 
is best understood by looking at the g evolution of the 
rapidities as shown on Fig. 2 for the three states in Fig. 
1. 

As g rises, each rapidity has a tendency to go down 
towards —00. Any block of unoccupied states stops 
up to Hi rapidities from doing so by keeping them fi- 
nite. If Pi rapidities are going down and they meet a 
block of Hi unoccupied states every rapidity will be kept 



S units of the interievel soaclng] 




FIG. 2: Evolution of the rapidities (real part) from to large 
g for the states presented in Fig. 1. 



8 



finite if Hi > Pi {Ai = rapidities will go through). 
On the other hand, whenever Hi < Pi, only Hi rapidi- 
ties can be kept finite and the remaining Ai = Pi — Hi 
will keep going down towards — oo. Starting from the 
highest block of rapidities (of size Pi) we therefore 
have Ai = Pi — Mm{Pi, Hi) which go through the 
Hi empty states below. These will be added to the 
following block of P2 rapidities giving P2 + Ai rapidi- 
ties which then meet a block of H2 unoccupied states. 
A2 = {P2 + Ai) - Min(P2 -I- Ai,H2) wih go through and 
continue their descent. Keeping this analysis going until 
we reach the last blocks gives out the result given above. 



2. Application to ground states and single excitation states 

At 5 = 0, for any of the canonical ground states con- 
taining Np rapidities in the Np lowest energy levels, we 
have a single block of unoccupied N — Np levels and a 
single block of Np occupied levels (see the left group in 
Fig. 3). Since no unoccupied levels are below this Np- 
rapidities block, they would all diverge in the g — ^ 00 
limit. 



0000000 
0000000 
0000000 
0000000 



Nr = Np 



Nr= Np-1 



Nr = Np-1 



FIG. 3: A ground state and the corresponding set of states 
which give a single finite rapidity at jr — >■ 00. 



because of Hi, and one of the P2 rapidities will remain 
finite as well, leading to a total of two finite rapidities. 

From very simple combinatorics we therefore find that 
we can build A'^— 1 single finite rapidity states (at g — 00) 
by deforming 5 = two-particle states. 

Moreover, at 5 cx) for a given number of rapidities, 
the total number of states with J = Np — 1 (single finite 
rapidity) is given by the total number of solutions to the 

single Bethe equation J2iLi x=b" = This equation 
also has A^ — 1 distinct solutions and therefore every state 
with a single finite rapidity at strong coupling stems from 
one of the two-particle states at 5 = 0. This was pointed 
out before in Refs. 30,31. 

From the previous subsections we concluded that for 
the form factors one can get contributions coming 
from the Np rapidities states with either A'p or Np — 1 of 
them diverging. The first case we now know corresponds 
to the ground state. Since we also showed that every 
state with one finite rapidity is generated by singly ex- 
cited states, it becomes clear that the two-particle states 
do give out every non-zero contributions in the g ^ 00 
limit. For S~ we showed that only the A'p — 1 ground 
state (all rapidities divergent) or the Np — 1 states with 
A'p — 2 diverging rapidities contribute. Once again this 
means that the intermediate sum can be limited to the 
Np — 1 ground state and single excitation (two-particle) 
states. One should also notice that, in this limit, any 
single excitation state which leads to two finite rapidi- 
ties will not contribute although at weaker coupling they 
could. 

Since the complete set of two particle states (plus the 
ground state) saturates the sum rules in both the .g — > 
and g — >■ cxD limit, it is reasonable to assume that they will 
also be largely dominant in the crossover regime. This 
fact will be explicitly proven numerically since even for 
the smallest systems, this subset represents, for any g, 
more than 99% of the weight. 



Focusing on the states built out of a single excitation 
above this ground state, we find 3 distinct cases. When 
one moves the top rapidity (from level Np) to a higher 
level (say level (3 > Np), the resulting block structure is 
(Pi = l){Hi= (3- Np){P2 =Np-l) (see middle group 
in Fig. 3) which leads in the strong coupling limit to the 
single Pi rapidity being kept finite while the remaining 
Np — 1 will diverge. In the second scenario, if one pro- 
motes any of the rapidities at level a < Np (excluding the 
topmost one) to the level A'p -I- 1 the resulting structure 
is {Pi=Np-a + l){Hi = 1)(P2 = a - 1) (see the right 
group in Fig. 3). Once again, this leads to one of the 
Pi rapidities staying finite due to the unoccupied block 
Hi ~ 1 while the other Np — 1 rapidities will diverge. 

The last possible case is a single excitation obtained 
by moving a rapidity from level a < Np into an empty 
level /3 > Np + 1. Doing this gives the following block 
structure (Pi = l){Hi = /3- Ap - 1)(P2 = Np-a){H2 = 
1)(P3 = a—1). In the end the Pi rapidity will stay finite 



V. DYNAMICAL CORRELATION FUNCTIONS 

We have presented all the ingredients to calculate the 
dynamical correlation functions: we need to solve the 
Richardson equations for each state calculate its 

energy E^, use the rapidities defining the solution to 
compute the determinant and calculate the form factors. 
However, while the formulas we obtained for the corre- 
lation functions are completely general and are valid for 
any choice of the Hamiltonian parameters Eq, and g, to 
obtain a physical result we still have to perform the sum 
over the states and this cannot be done analytically. Thus 
we need to make a choice of the model to study. As we al- 
ready mentioned, we only consider the most-studied case 
in the condensed matter literature, which consists of A^ 
equidistant levels at half-filling, i.e. A^ = 2Np. We define 
the levels as 

ea=a with a^l...N , (42) 
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FIG. 4: Diagonal correlation function G'^^{lo) obtained by 
smoothing the energy (5-function with a Gaussian of width 
WQ = 0.1 for different values of g and at fixed number of 
pairs Np = 32 [N = 64). We show C{g) + Gi^{uj){10^)o^ 

_9 1 

0.05 

with the offset C{g) = ^ 20(10^)" (Effectively this makes 

n = 

scaled logscale plots with an offset). 



i.e. we measure the energy scale in terms of the inter- 
level spacing and we fix the Debye frequency (the largest 
energy level) to N. 



A. Diagonal Sz correlator 

We start our analysis with the diagonal Sz correlator 
that, in frequency space, reads 

N 



G'^ (a.) - W l(Wl^n|g^)P ... 



Ev + Eos) , 
(43) 

with |{t;}) having N/2 rapidities. In particles language 
this is a density-density correlator. At any finite N, this 
is a sum of 6 peaks each at the energy of the excited 
state [{t'}), and each weighted with the corresponding 
form factor. 

In the thermodynamic limit where N ^ oo while the 
inter-level spacing d — > (keeping a finite bandwidth 
for the single particle excitations), we can think of any 
5 > as being already the strong coupHng case, i.e.: the 
BCS mean-field treatment becomes exact. In such a case 
correlation functions should be described by the 5 — ?> oo 
limit where only the single finite rapidity band would con- 
tribute. However, we are interested here in mesoscopic 
effects at finite that are encoded in the quantization 
of the energy levels and in the potential presence of non- 
trivial contribution from other excited states. 

In any real physical system, the 5-peaks are smoothed 
by different effects such as temperature broadening or, 
even in the case of very small T, experimental resolution. 
All these effects are expected to broaden the (5-peaks in 
an approximately Gaussian fashion. For this reason, in 
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FIG. 5: Diagonal correlation function 0^^(10) obtained by 
smoothing the energy 5-function with a Gaussian of width 
wg = 0.1 for different values of g and at fixed number of pairs 
Np = 32. We present the first excited subband contribution 
plotted with the frequency rescaled by the grand-canonical 
BCS gap, bringing the bottom of the band at cj/Ascs ~ 1. 



Fig. 4 we plot a typical example of such correlator for 
different values of g, at fixed number of pairs Np = 32, 
with the 6 functions broadened to Gaussians of width 
Wg ~ 0.1. The necessity to use logscale plots comes 
from the fact that typically the second energy subband, 
clearly separated by a gap for g > 0.45 has contributions 
which are orders of magnitude below the contributions 
from the first excited subband. In the limit 5 — )■ 00, the 
width of these two subbands would go to zero while the 
gap separating the first subband from the ground state 
and the second band from the first, would tend to the 
grand-canonical BGS gap: 



N 



2sinh(l/2g) 



(44) 



The quantization of the energy levels is evident, espe- 
cially for small g. For larger values of g, a peak seems 
to develop at the BCS gap even for these relatively small 
values of N. The contributions from the first excited sub- 
band are made more apparent in Fig. 5 which shows only 
this subband. We can understand this peak developing 
at the bottom of the band as being mostly due to the 
growing density of states at this energy. It is not due to 
a single large contribution from a given state, but sim- 
ply to the finite width of the gaussian peaks making the 
small energy differences unresolvable. In fact, as shown 
explicitly in appendix B, in the strong coupling limit all 
states in the subband have identical form factors and 
therefore equal contributions to the correlation function. 
The coalescence of the energies happening faster (in 17) 
at the bottom of the band lead to this impression of a 
developing peak. 

Very similar plots can be found frequently in the ex- 
perimental literature (see e.g. the review [3]) for the I-V 
characteristic of superconducting nano-grains. We are 
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now in a position to understand these results quantita- 
tively. 

However, without a realistic idea of the kind and am- 
plitude of broadening these plots are still only indicative. 
Moreover, it is difficult to correctly visualize these func- 
tions and the relative importance of the various contribu- 
tions. A more precise information about the mesoscopic 
effects is encoded in the integrated correlation function 
obtained from Gf^ integrating it up to a given a; 



duj'Gtiu:'), 



(45) 



i.e. the sum of the form factors of states with energy 
smaller than lo. The integrated correlation functions can 
be plotted without any smoothing. For several values of 
g, we report them in Fig. 6 for fixed N = 16. Any step 
corresponds to a different eigenvalue (most of them have 
a two-fold degeneracy) of the Hamiltonian. Notice that 
for large uj the integrated correlation function gives the 
value of the sum rule being the sum of all form factors, 
i.e. I^iu; = oo) = {(S^,,)^) (cf. Eq. (27)). ^ 

To understand how the various excitations combine to 
give this correlation function, in Fig. 7 we report the 
contribution of the ground state and of the two-particle 
states to the sum rule. The ground state (inset) ac- 
counts for the full correlation function in the absence of 
interaction, (as discussed previously) but its contribution 
quickly decays to zero with increasing g, as consequence 
of a complete re-organization of the ground state struc- 
ture. In the main plot of Fig. 7, we show the sum of 
the ground state plus all two-particle states. It is evident 
that in all considered cases, these states give basically 
the full correlation function, and at most about 1% is 
left to the other states. In particular for large g, the two- 
particle states account for the full correlation function 
(as shown in the previous sections). Notice also that with 




FIG. 6: Integrated diagonal correlation function I^zi'^) for 
= 16 and varying g. The inset is a zoom of the upper 
left corner, where the quantization of the energy level is more 
evident. 




1 1.5 2 

g [in units oF th& interl&wsl spacing) 



FIG. 7: Sum-rule for the diagonal G'^^{lo) correlation func- 
tion. Inset: The ground state contribution. Main Plot: 
Ground state plus all two-particle states. 



increasing N , the missing contribution becomes smaller. 
This property simplifies enormously the computation of 
the correlation functions. In fact, to have an effective 
description of the correlation functions, we do not need 
to sum over the total (^2) states, but only over the 
N'^/4: two-particles states. For this reason, in the vari- 
ous figures, all correlation functions have been calculated 
by only considering two-particle states (for lower N, we 
checked by full sums that this 'approximation' does not 
introduce any visible modification.) 

We can now come back to the analysis of the corre- 
lation function I^zi^) itself. In Fig. 8 we plot three 
different curves at fixed gN in such a way that the BCS 
gap is constant in the large g approximation (see Eq. (44) 
giving the BCS gap in the grand-canonical ensemble). It 
is evident that for the smallest value of = 16, which 
corresponds to the largest value oi g — 2.8, /^^(w) has al- 
most its asymptotic expression (see Appendix B) that is 
a step function at the BCS gap (the integral of S{Ld — A) 
contribution for large g) . Although the first energy band 
still has a finite width (which would go down to zero as 




FIG. 8: Integrated diagonal correlation function Izzi^^) for 
three values of A'^ and g, while keeping constant Ng. 
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FIG. 9: Integrated diagonal correlation function /^^(a;) for 
A'' = 64 and various g. The horizontal axis has been rescaled 
by the BCS gap. 



according to equation B13) this width is akeady sig- 
nificantly smaller than the BCS gap (which scales as g). 
Oppositely for TV = 64 and g = 0.7, despite N being 
larger, the small value of g and the still significant rela- 
tive width of the sub-band leads to a nice quantization 
of the energy spectrum, resulting in a staircase function 
with unequal steps. In the inset of Fig. 6, the zoom of 
the lower energy sector for small 5 at iV = 16 is reported, 
showing the formation of intermediates steps of unequal 
contributions with increasing g. 

Finally in Fig. 9, we report /^^(w) at fixed = 64 
and varying g, rescaling the horizontal axis by the grand- 
canonical BCS gap (eq. 44). 

For large g, Aqc corresponds exactly to the energy 
of the first excited state (and therefore sub-band) and 
all the curves start rising from 1. However for smaller 
values of g (in particular for g — 0.5), visible finite size 
effects are present. As already stressed above, the stair- 
case structure for small smoothly connects to a step 
function for large enough g. 



B. Diagonal {S^S ) correlator. 

This correlation function (as well as the corresponding 
off-diagonal one) is directly connected with the annihi- 
lation and creation of a Cooper pair. We could imagine 
for example an experiment in which a superconducting 
grain was contacted on the left and right by two separate 
bulk superconducting leads, each having their respective 
order parameters, a problem which was studied for exam- 
ple in [33,34]. The Josephson current in such a system is 
given by perturbation theory in the coupling between the 
leads and the grain, and its calculation involves comput- 
ing such correlators as the one above. We will treat this 
problem in more detail in a separate publication, concen- 
trating for the moment on the results for the correlators 
themselves. 




0.5 1 1.5 2 2.5 3 

g [in units oF th& interlffwsl spacing) 



FIG. 10; Sum-rule for the diagonal G+_(a;) correlation func- 
tion. Inset: The ground state contribution. Main plot: 
Ground state plus all two-particle states. 

We start by showing the sum rule (c.f. Eq. (28)) ob- 
tained from the ground state and all the two-particle 
states as function of g at fixed N in Fig. 10. For the 
ground state, we have an opposite behavior compared to 
the correlation: the value for small g is low and it 
increases by increasing g, but never saturating it com- 
pletely. The missing contribution is again mostly in the 
two-particle states, as evident from Fig. 10. For 5 = 
and g — > oo, it has been proved analytically in the pre- 
vious section that this class of states gives a perfect sat- 
uration of the sum rule and hence the full correlation 
function. The figure shows that even for intermediate g, 
a very accurate calculation comes only from two-particle 
states, and so in the following we will ignore all the other 
too small contributions. Notice that the contribution of 
the leading excitations is even more relevant than for Gf^ , 
having always saturation above 99% and that again for 
large enough N it increases while increasing N. 

The integrated correlations 

I^_{oj) = / du'G'l_{oj'), (46) 
Jo 

are reported in Fig. 11 for N — 64. In the left panel we 
report the small g < I. At very low g, the BCS gap is 
not yet formed and the correlation function has a stair- 
case structure with almost equal steps, as a consequence 
of the almost perfect equispacing of the energy level. In- 
creasing g, the levels become incommensurate and the 
correlation function acquires a structure that reflects the 
formation of the BCS gap. In the right panel we show 
the same correlation function plotted in terms of uj/Agc 
with the grand-canonical gap given by Eq. (44). For 
large values of g, the 'band structure' of the energy lev- 
els is evident in the formation of two main-steps at Aqq 
and at 2Agc- In the same panel we also show some 
intermediate and small values of g to show the forma- 
tion of this two-step structure from the sharpening of the 
smaller steps. This is the main signature of mesoscopic 
effects in the correlation functions. Whereas for the 
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correlations the second excited band's contribution were 
very rapidly supressed, in the case at hand, it maintains 
a very important contribution in the full regime stud- 
ied in this paper. Although the selection rules point out 
that it no longer contributes at g — ?> cxd it appears that 
the suppression of the form factors happens at much 
stronger coupling than it does for . This seems to sug- 
gest that corrections to the mean field BCS results (more 
or less equivalent to the g ^ oo limit) could show up in 
a wide regime when looking at the previously mentioned 
Josepshon current experiments. 



C. Off-diagonal correlators. 

The off-diagonals correlators have properties very sim- 
ilar to the diagonal ones. For this reason we only shortly 
present them, without an extensive discussion. In Fig. 
12, we report the contribution to the sum-rule (c.f. Eq. 
29) from the ground state and from all two-particle 
states. The ground state by itself saturates the sum rule 
for large enough g, while for small values the contribu- 
tions of the two-particle states are essential. Notice that 
the saturation of the sum-rule is always above 99.9%, 
even more than in the diagonal case. Since the ground 
state is the strongly dominant contribution to the corre- 
lation, G°^_ will be a (5 peak at the ground state value. 
Thus this is the less sensitive function to reveal meso- 
scopic effects. 

In Fig. 13, we report the integrated correlation func- 
tion 

/f_H= / doj'GlU^'), (47) 
Jo 
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FIG. 11: Integrated correlation function as function 

of uj for several couplings g and fixed = 64. Left: small 
values of (? < 1. Right: as function of cj/Agc to show the 
formation of two steps in the limit of large g. 
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FIG. 12: Sum-rule for the global G°^_ (cj) correlation function. 
Inset: The ground state contribution. Main plot: Ground 
state plus all two-particle states. 



for several g and N = 64. To make the plot visible we di- 
vide by the off-diagonal order parameter od as obtained 
by the Hellmann-Feynman theorem^^. In the main plot, 
small g < 1 are reported. As before, for small enough 
g we have a regular staircase behavior, that however be- 
comes quickly fiat because the ground state contribution 
to the correlation function is too large. For large values 
of 2.8 < .g < 3 are reported in the inset of the figure 
and are zoomed very close to 1. The contribution from 
two bands at I^gc and 2/S.qc is evident, but it is visi- 
ble only because we have an exact solution at hand. It 
would have been very hard, if not impossible, to see such 
a small effect in any numerical approach and in a real 
experiment. 
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FIG. 13: Integrated diagonal correlation function I'^_{lo). 
Main plot: Small values oi g < 1 at step of 0.05 going 
from the staircase to the flat regime. Inset: Large values 
oi g = 2.8, 2.85, 2.9, 2.95 showing the very small two steps 
structure. 
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FIG. 14: Real-time correlation function with varying g. Top: A\G'lz(t)\/N; at small to intermediate coupling (left panel); as a 
function of rescaled time in the strong coupling regime (right panel) . 

Bottom: 2\G'!^_{t)\/N\ at small to intermediate coupling (left panel); at strong coupling (right panel). These functions have 
been shifted by an integer to make it easier to read, but they remain bounded between and 1. 



D. Real time correlations. 

It is also interesting to look at these correlation func- 
tions in real time. The presence of incommensurate en- 
ergy levels gives indeed a quite complicated structure 
that greatly simplifies in the best known large and small 
g behavior. We show in Fig. 14 the typical evolution of 
the two diagonal correlation functions. It is evident that 
for small 5, the functions are almost periodic as a trivial 
consequence of the almost commensurate levels (this is 
strictly true only at g = 0, but at g = 0.1 several peri- 
ods should pass before the non-exact commensurability is 
manifest). Increasing the behavior becomes irregular 
and for very large g, new structure emerges due to the 
formation of a well defined energy subband structure. 

As can be seen on the top right panel, for large enough 
(7, we find a peculiar scaling behavior for the norm of the 
Gf^ correlation function. Indeed, dividing time by 5', 
the dimensionless coupling constant {g' = g/d with d the 
interlevel spacing) we find an almost perfect agreement 
between the curves for various strong coupling cases. 
The specific details allowing a clear understanding of this 
regime are presented in Appendix B where we compute 
the \/g' expansion of both the energies and the form fac- 
tors. One should know, however, that looking indepen- 
dently at the real and imaginary parts of this correlation 
function they both still show a rapidly oscillating compo- 



nent associated with the gap frequency. It is only when 
looking at the norm of the correlation that this surprising 
scaling can be found. 

The same is not true for the Gi^_ correlator shown on 
the bottom two panel of figure 14. In this case, even 
when looking exclusively at the norm of the correlation a 
contribution oscillating at the very large energy gap fre- 
quency remains present. Due to the time scales plotted, 
this rapid oscillation can barely be resolved on the figure 
for coupling strength larger than 0.40 but a rescaling of 
the plots clearly shows that they remain present. 



VI. CONCLUSIONS 

We have studied the dynamical correlation functions 
of the reduced BCS model in the canonical ensemble by 
means of Algebraic Bethe Ansatz techniques. We pre- 
sented analytic selection rules in the weak and large cou- 
pling regimes. For finite values of the coupling and for 
finite number of particles, the correlation functions are 
calculated numerically by summing over the form fac- 
tors of the relevant states. We showed that two-particles 
states always gives saturation of the sum rules that is 
above 99% for all g and N considered, in stark contrast 
with other integrable models like spin-chains and one- 
dimensional gases. We presented and discussed exten- 
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sivcly the crossover from small to large g, where quantum 
fluctuations have a dominant role and lead to a differ- 
ence between the canonical and grand-canonical ensem- 
ble, and where a perturbative treatment would provide 
incorrect results. 

Additionally, we showed unexpected behavior of 
correlations at half-filling. The properties of these cor- 
relations give rise to a scaling law which is the exact 
opposite of what one would expect if the properties were 
dominated by the energy scale associated to the super- 
conducting gap. 

Wc remind the reader that while wc only studied the 
case of N non-degenerate equidistant energy levels at 
half-filling (which is the most interesting model from the 
condensed matter point of view^'^"''), this is in no way a 
strict limitation of our approach, which can be in princi- 
ple applied to any choice of initial energy level distribu- 
tion. The description of pairing in nuclei is for example 
a situation where other choices of the parameters are 
more natural*'"''^"'''**. These could be treated by a simple 
adaptation of our results. 

The correlation functions we have provided have ap- 
plicability in the phenomenology of many experimen- 
tal situations. In further work, we will apply the cur- 
rent results to transport phenomena through metallic 
nanograins coupled to normal and/or superconducting 
leads. 
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Appendix A: Solving Richardson equations 

As explained in our previous publications^^'^^'^'^ (see 
also [30,39-42]), we find the solutions starting from the 
trivial g = solutions and slowly raising the value of 
the coupling constant. With a linear regression using the 
previously found solutions, one gets an appropriate guess 
for the new solution at g + Sg. Using a simple Newton 
algorithm for solving coupled non-linear algebraic equa- 
tions is then sufficient, provided one makes the necessary 
change of variables. 

Indeed, at any value of the coupling, rapidities arc ci- 
ther real or form complex conjugate pairs (CCP) and the 
pairing of two rapidities only occurs at bifurcation points 
g* at which they arc both exactly worth ec, i.e. one of 
the single particle energy levels. It is also possible, as g is 
increased, that two paired rapidities split apart becom- 
ing both real again. At a critical g* where a pair splits 
or forms, the derivatives ^p- are not defined making the 



computation of the necessary Jacobian impossible. This 
problem is easily circumvented by, in the vicinity of crit- 
ical point at which Wi = Wj = e^, making the following 
change of variables 

A+ = Wi + Wj 

A_ = {wi-Wjf. (Al) 

On both sides of g* , those are two real variables and 
they have well defined derivatives even at the bifurcation 
point. The fact that g is slowly increased allows us to 
figure out beforehand whether given rapidities are about 
to form (or break) complex conjugate pairs. Naturally, 
it makes the numerical procedure more tedious than it 
would be if one was able to guess correctly the struc- 
ture at the precise value of g in which we arc interested. 
However the lack of known analytical results about the 
solutions to these precise Bethe equations forces us to use 
this scanning procedure. Fortunately, for the dynamical 
correlations we only need a very restricted set of states 
in order to get a very accurate description. Since single 
solutions are addressed one by one independently of the 
dimension of the full Hilbert space, the problem remains 
numerically tractable for fairly large system sizes which 
matrix diagonalization could not tackle. 

In figure 15, we show the real part of the Np = 16, N = 
32 rapidities as a function of g for typical single excitation 
states. As in the rest of the paper the lowest single parti- 
cle energy level is chosen to be ei = 1 and the subsequent 
are at — a. 

For the ground state, rapidities form complex conju- 
gate pairs in a very simple fashion, the top one pairing 
with the next one below and so on. For an odd number 
of rapidities, the lowest one would remain real and the 
other Np — 1 form CCPs in this way. 

In the top right pannel, we show a state obtained by 
promoting to level (3 > Np + 1 the highest rapidity (level 
index a = Np) from the ground state. We know from the 
algorithm presented in section IV C that this state will 
have a single finite rapiditiy at strong coupling. More- 
over, one sees that the promoted rapidity will simply stay 
real and remain between and e^_i. Since there is 
no other rapidity close by, it cannot form a pair to go 
through the energy level e/3_i. Indeed, the structure of 
the Richardson equations prevents a single rapidity to be 
equal to an energy level eg since the diverging term 

has to be cancelled by a diverging term ~^ . The re- 
maining rapidities will simply form CCPs as they would 
for a A^p — 1 pairs ground state. 

The two figures in the middle are built by promoting 
a rapidity with index a < Np to level /3 = Np + 1. In 
both cases, we find only a single finite rapidity at strong 
coupling but we still have two different scenarios. The top 
contiguous block has an even (odd) number of rapidity on 
the left (right) pannel. For an odd number of rapidities in 
the top block, the lowest one from that block will simply 
stay real and between ea+i ^nd £„. In the even case, 
the lowest two rapidities from the top block will actually 
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g (in units of the interlevel spacing) g (in units of tlie interlevei spacing) 




g (in units nfttie interlGvel spacing) g (in units cfttie intarlevei spacing) 



FIG. 15: Real part of rapidities for: ground state (upper 
left), promoted top (q — Np) rapidity (upper right), rapidity 
promoted at p — Np + 1 right above the Fermi level (middle 
left), rapidity promoted a.t P — Np + 1 right above the Fermi 
level (middle right), non-top rapidity promoted above A'j, + 1 
level (lower left), non-top rapidity promoted above Np + 1 
level (lower right) 



Appendix B: Strong coupling scaling 

Hand-waving arguments tend to lead to the idea that 
the excitation gap should be the dominant energy scale 
in this system at least for strong coupling where one ex- 
pects the BCS description to be more or less adequate. 
We would therefore expect correlations in time to show 
strong oscillations at a frequency given by the gap (i.e. 
oc g in the strong coupling limit). Surprisingly, the S^S^ 
correlations shown in figure 14 actually show the exact 
opposite behavior. The magnitude of correlations do 
show scalable behavior but they happened to be slowed 
down by an increasing g and therefore an increasing gap. 

This appendix aims at understanding this peculiar 
scaling property. In order to do so, we can rely on the 
strong coupling expansion of the Richardson equations. 
The following analysis will give us 'semi-analytical' (i.e. 
getting the numerical values still needs numerical work) 
expressions for both the energies and the form factors 
needed to understand this correlator. 

We assume a convergent expansion exists around the 
g ^ oo point for the values of the rapidities themselves 
(Fig. 15 shows this assumption to be correct at strong 
enough g). Keeping only the first relevant corrections, we 
can write down, for the 3 types of states we are interested 
in: 

the ground state 

Wj = Cf^g + Aj +Bj- \J j = l...Np (Bl) 
9 

the single finite rapidity states 



form a pair at the energy level a -I- 1, then split apart into 
two real rapidities at level a. This splitting will make 
one rapidity be above where it will stay until g ^ oo 
while the other one goes below and will later form a 
pair with the top rapidity from the block underneath. 

Finally, the lowest pannels which lead to two finite 
rapidities at strong coupling behave in a similar way for 
the lowest Np — 1 rapidities, whereas the promoted one 
stays finite in the same way it did in the top right pannel. 



At (7 — >■ oo the state with one single finite rapidity in 
between e^_i and e^, for A^p < f3 < N is the deformed 
version of the 5 = state built by exciting the ground 
state's top rapidity Wi = epf^ to level ep. On the other 
hand if the finite rapidity is between 1 < a < Np and 
a + 1 the state is the deformed version oi g — state ob- 
tained by promoting the ground state's rapidity Wi = 
to level e^Vp+i- The approximative location (between 
given energy states) of these strong coupling finite ra- 
pidities was already know and used in^^ to compute non- 
equilibrium dynamics in the related central-spin model. 
However we here show how they each correspond to a 
known g = state deformed by interactions. 



1 

9 



C. 



1 V 3 



\...Np - 1 (B2) 



9 



and states with two finite rapidities 



V2.k = »/Vc +<52,fc- 



3 1 

J J g 



V j = l...Np - 2. (B3) 



In the last two cases we have respectively either — 1 
or {NpY — NpN — N — 1 possible values of k each asso- 
ciated with a different possible solution of the Richard- 
son equations (different possible values of the i] rapidities 
which remain finite). In all three cases, eq. 31 gives the 
values of the C constants which define the diverging ra- 
pidities in the g — >■ oo limit which are therefore indepen- 
dent of the choice of finite rapidities values for any given 
k. One should understand that g here is considered to 
be the dimensionless quantity ^, d being the interlevel 
spacing. 
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1. Energies 

The respective energies of these states, obtained by 
summing the Np rapidities are therefore simply given by: 



(B4) 

















■3 = 1 


9 + 


.7 = 1 


+ 


E^^ 



These equations arc defined for any of the Np — 1 values 
of index j. Summing up the Np — 1 equations B9 (divided 
by Cf""^), we find 



~{Np-l) = iN-2)J2-^, (Bll) 
while summing up the equations BIO gives us: 



Ek = 







E c?-' 


9 + 







vr+ E ^: 



+ 



Np-l 



(B5) 
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(Np - 1) 
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(B12) 
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Np-2 
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(B6) 



The expansion of the Richardson equations around the 
strong coupling solutions with a single finite rapidity 
gives us: 



— 



N 



1 



ATp-l 



a— 1 '^j J T^J 



-2^B7) 



Afp-l 



- 2) - 2 E ^^p-i ^iVp-i 



E' 



iVp-1 /^Wp-l , 

'r -^3 



1 



(B8) 



Which, order by order gives: 



Using this last expression, the energies of every single 
finite rapidity state can be written including the lowest 



correction in - as: 



Ek = 



N„-l 



E^; 

Np-l 



N„-l 



This shows that the half-filled case leads to the com- 
plete energy collapse of the first excited band as proven 
before in [27], i.e. lim Ek — Ek' = 0. The zero band- 



g->-oo 



width obtained in this specific case will be shown to be 
one of the central elements in the scaling properties of 
the 5^ operators dynamical correlations. 



2. Eigenstates 

In order to establish a similar expansion for the states 
themselves, one can simply use their representation as a 
Bethe state (eq. 11) and expand the C operators used to 
construct them. Keeping terms up to order ^ for both 

divergent {w = Cg + A + ^) a,nd finite [t] = r]°° + ^) 



rapidities we have: 



_^Wp-l 



-A? = 



■ N 

E' 



„JVp 



-^V°k 



C. 



JV„-1 



■'Vp-1 l(J^p-'^^k _ j^k (J^P^^) 
+2 E !' .. f i (BIO) 
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Civ) 
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JU -Ca 
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gc'''"*^gC 
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E^^ 
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c g 



(B14) 
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By getting rid of the ^ prefactors which, for physical 
quantities, will always be cancelled by equivalent factors 
in the norms, we can therefore write: 



N 



\GS) « \GS^'') + -J2GaS+\GS^^-^) (B15) 



N 



1%) « 

a=l 
1 ^ 

+- J2 Gl^S+S+\GS''^-') (B16) 



N 



1 ^ 

+- E G^^p ,^S^S^S^ IGS^" ^) 

(B17) 



a,jS,7=l 



where the states \GS^) = (S'jtt)^ |0). The following set 
of definitions was also used: 



Ga 
pk 



^ A- 

.7=1 '^3 
1 



.7 = 1 '-".7 



^Ar„-i 



4/c 

Wp-i ^ 

.j'=l 
1 1 



1 



'52,fc 1 



done previously (cq. 39). Although the first order cor- 
rection can also be obtained in a similar fashion, we will 
not explicitly need the coefficients of the expansion and 
therefore simply write: 



{GS\Sl\GS) 



N -l\ fN -I 



Np-l 



+ A-a, 

9 



(B19) 



with Aa an unspecified (although obtainable) constant. 
Identically, for single finite rapidity states we can write: 



N-2 • 



a)® Il"^"^ •■•taNp-2|^ 

{ai, ... ajVp-2} 

\Np-lJ 

"o It/34 



{ai, ... awp-i} 



/3#Q 



(jVp-l) 
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iV-2\ /iV-2 



Np-2J \Np - 1 



N-2\ /N-2 



2 niv„-i 



N-2\ /N-2 
N„-lj \Np-2)^\N,,-l 



N -I 



(B20) 



The orthogonality of these states with the ground state 
also allows us to write: 



1 '-'.1' 



^ -J- J 
■iVp-2 ^ 

.7' = 1 '^.7' 



(B18) 



(G5 iT^fe) 



i.Np-i; 

It/S)® E |{tai - 



3. Form factors 



N-1 
Np-1 



{ai, ... awp-i} 
0, 



(B21) 



At g ^ 00 it is straightforward to compute value of and therefore, adding the next order term through an 
the various forms factors. For the ground state it was unspecified constant: 
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pk 



N-1 
Np-1 



N-2' 



+ - 



{Np - iy.{N - Np - 1) 



N-2 
Np-1 



+ (B22) 



(the ones involving the single finite rapidity states) arc 
actually all equal. At the next leading order the only 
contributions also come from the same reduced set of 
states. 

Using equation 43 we can write the correlation function 
by summing over the A''— 1 possible values of r}k- At order 
we have 



It was also proven in section IV B that at g — > cxd wc 



have {GS\ \m,k,m,k) 



g->-oo 



and we therefore have: 



{GS\S^ \r]i,k,r]2,k) 



(B23) 



Finally one can similarly compute the squared norms 
of the ground state and the single rapidity states. 
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iVklVk),^^ = Y^F^iF^r {GS'''^-'\S^S+\GS''''-') 
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Specializing to the half-filled case, we find: 
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At order in ^ wc therefore find that the only A^ — 
1 non-zero contributions coming from the form factors 



N-l 



GUu;)^J2^(''-^k + EGs) 



k=l 



N 

4(A- 1) 



N 



2Re [f:;B^\ {N/2)\{N/2)\ 



+E 



whose Fourier transform gives us 
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Looking exclusively at the magnitude of the correla- 
tions and therefore at phase independent properties of 
this correlator we have 
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Since equation B13 tells us that 



Ek. — Ew = 
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we can write 
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This shows that at strong enough coupling the domi- 
nant term 

is purely a function of |. This fact is only true at half 
filling where the ground state average of any S^^ is zero. 
Moreover, half filling also allows fulfillment of the second 
necessary condition, the vanishing of width of the first 
excited band. The first effect makes the energy scale as- 
sociated to the BCS gap irrelevant since only the first 
band of excited states is contributing to the correlations. 
The vanishing width of this band is then responsible for 
the 'inverted scaling' by making the relevant energy dif- 
ferences smaller as g increases. 



This very peculiar half-filling scaling property makes 
it possible to slow down specific dynamical processes in 
this system by making the interaction stronger. Since 
they are only related to the spectrum and the condition 
{GS\ \GS) = 0, this scaling law would hold at half- 
filling for any possible correlations of operators; be 
they local, global, intra or inter-level correlations. The 
magnitude of any one of the possible correlations would 
still follow a similarly scalable time evolution. 

Although clearly valid in the strong g limit where ^ < < 
1, we also find through the numerical work carried out 
in this paper that this scaling behavior extends to a very 
broad range of coupling constants. For g > l.bd with 
d the inter-level spacing, we find that ^ corrections are 
already strongly suppressed and the scaling behavior is 
therefore already apparent. 
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